#ifndef xpack_math_cos
#define xpack_math_cos
#pragma push_macro("xuser")
#undef  xuser
#define xuser mixc::math_cos::inc
#include"define/base_type.hpp"
#include"math/abs.hpp"
#include"macro/xalign.hpp"
#include"macro/xexport.hpp"
#include"math/literal.hpp"
#include"math/private/parallel_lut.hpp"
#include"utils/array.hpp"
#pragma pop_macro("xuser")

namespace mixc::math_cos{
    /*
    void gen(){
        f80 x = 1;
        for(uxx i = 0, j = 0; i < 10; i++){
            x *= ++j;
            x *= ++j;
            printf("type_t(%+.18Le), \n", f80(i & 1 ? +1.0 : -1.0) / x);
        }
    }
    */

    template<class type_t, uxx multiple>
    inline xalign(32) auto lut = inc::parallel_lut<multiple>(
        type_t(-5.000000000000000000e-01), 
        type_t(+4.166666666666666667e-02), 
        type_t(-1.388888888888888889e-03), 
        type_t(+2.480158730158730159e-05), 
        type_t(-2.755731922398589065e-07), 
        type_t(+2.087675698786809898e-09), 
        type_t(-1.147074559772972471e-11), 
        type_t(+4.779477332387385298e-14), 
        type_t(-1.561920696858622646e-16), 
        type_t(+4.110317623312164858e-19)
    );

    namespace inc{
        using namespace literal;
    }
}

namespace mixc::math_cos::origin{
    // 注意：
    // 此 unsafe 函数不带定义域检查，需要使用者保证 x 的范围属于 [-pi, +pi] 的闭区间
    template<class type_t, class ... args_t>
    inline auto cos_unsafe(type_t const & x0, args_t const & ... xn) {
        constexpr auto multiple = 1 + sizeof...(args_t);
        constexpr auto step     = sizeof(type_t) == sizeof(f32) ? 5 : 9;

        auto && x2      = inc::make_array((type_t &)x0, (args_t &)xn...);
        auto && xm      = inc::array<type_t, multiple>();
        auto && rx      = inc::array<type_t, multiple>();

        for(uxx i = 0; i < multiple; i++){
            x2[i]      *= x2[i];
            xm[i]       = x2[i];
            rx[i]       = 1;
        }

        for(uxx i = 0; i < step; i++){
            for(uxx j = 0; j < multiple; j++){
                rx[j]  += xm[j] * lut<type_t, multiple>[i * multiple + j];
                xm[j]  *= x2[j];
            }
        }

        if constexpr (multiple == 1){
            return rx[0];
        }
        else{
            return rx;
        }
    }

    template<class type_t, class ... args_t>
    inline auto cos(type_t const & x0, args_t const & ... xn) {
        // note:
        // 直接用公式 cos(x) = sin(x + pi/2) 套壳在 x + pi/2 >= pi/2 时，
        // 会有 sin normalize 转换，导致只有一半性能，所以这里特化 cos 实现
        auto normalize = [](type_t x){
            constexpr auto rpi  = type_t(1.0) / (inc::pi);

            if (x = inc::abs(x); x <= inc::pi){
                return x;
            }

            auto q              = i64(x * rpi);                 // 整数倍
            auto r              = type_t(x - q * inc::pi);      // 余数
            return r;
        };
        return cos_unsafe(normalize(x0), normalize(xn)...);
    }

    inline f32 cos_unsafe(f32 x) {
        return cos_unsafe<f32>(x);
    }

    inline f64 cos_unsafe(f64 x) {
        return cos_unsafe<f64>(x);
    }

    inline f32 cos(f32 x) {
        return cos<f32>(x);
    }

    inline f64 cos(f64 x) {
        return cos<f64>(x);
    }
}

#endif

xexport_space(mixc::math_cos::origin)